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A general method for computing kinetic isotope effects is described. The method uses 
the quantum-instanton approximation and is based on the thermodynamic integration 
with respect to the mass of the isotopes and on the path-integral Monte-Carlo evalua- 
tion of relevant thermodynamic quantities. The central ingredients of the method are 
the Monte-Carlo estimators for the logarithmic derivatives of the partition function and 
the delta-delta correlation function. Several alternative estimators for these quantities 
are described here and their merits are compared on the benchmark hydrogen-exchange 
reaction, H+H2 — >H2+H on the Truhlar-Kuppermann potential energy surface. Finally, 
_C ' a qualitative discussion of issues arising in many-dimensional systems is provided. 
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C"| ■ 1 Introduction 

Measurement and theoretical predictions of kinetic isotope effects belong among 
,__! ' the main tools of chemical kinetics. Kinetic isotope effect (KIE) is defined as the 

^ . ratio JtA/ks of rate constants for two isotopomers A and B. Isotopomers A and 

B are two chemical species differing only by replacing a group of atoms in chemi- 
cal species A by their isotopes in species B. Recently, observation of anomalously 
large KIEs has helped prove importance of quantum effects in enzymatic reactions 
f— ^ . at physiological (i.e. surprisingly high) temperatures [1]. This and similar results 

\Q ' have changed our understanding of enzymatic catalysis and spurred an active ex- 

^^ . perimental and theoretical research in the last several years. 

r/3 ' Since the early days of chemical kinetics, KIEs have been predominantly de- 

scribed from the perspective of the transition-state theory (TST) [2,3]. This theory 
is intrinsically classical, although various quantum "corrections" have been incor- 
porated in it over time. These include corrections due to the zero-point-energy 
effects, high-temperature Wigner tunneling correction [2,3], and various semiclas- 
sical approximations for treating the tunneling at low temperatures [4-6]. On the 
other end of the spectrum are exact quantum-mechanical methods for computing 
\/ , rate constants and KIEs [7], but in general these are not feasible for systems with 

many degrees of freedom. One therefore resorts to various approximations that 
make a computation practicable but are less severe than the TST. Among these 
belongs a variety of quantum transition-state theories [8-10], the most recent of 
which is the quantum instanton (QI) approximation [11], motivated by an earlier 
semiclassical instanton model [4]. In this contribution, we describe a method [12], 
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based on the QI approximation, for computing KIEs directly, rather than via com- 
puting the rate constants for the two isotopomers first. Because of the ultimate goal 
of applying a similar method to enzymatic reactions, the method is implemented 
using a general path-integral approach that scales favorably with the number of 
degrees of freedom. Several alternative estimators for relevant quantities have been 
developed [12, 13] and their relative merits are compared in this contribution on 
the benchmark hydrogen atom-diatom exchange reaction. 

2 Quantum instanton approximation for the kinetic isotope effects 

The quantum-instanton approximation for the rate constant was introduced in 
Ref. [11]. A simpler alternative derivation [14] described in detail in Ref. [12], starts 
with the Millcr-Schwartz-Tromp formula [15] for the thermal rate constant k, 

/■OO 

kQ r = dtC s (t). (1) 

Jo 

Here Q r is the reactant partition function (per unit volume for bimolccular reac- 
tions) and Cff (t) is the symmetrized flux-flux correlation function, 

Cff(t) - tr (c-^/^c-^/V^Ffce- 1 ^) (2) 

with Hamiltonian operator H and flux operator F^ denoting the flux through the 
dividing surface 7 = a, 6. Quantum instanton expression follows by multiplying and 
dividing the integrand of Eq.(l) by the "delta-delta" correlation function C^it) 
defined below in Eq. (5), assuming that Cg(t)/Cdd(t) varies slowly compared with 
C'dd(t), and applying the steepest descent approximation to the resulting integral. 
Assuming further that the stationary-phase point is at t — 0, we obtain the QI 
thermal rate constant, 

^-i c »<°)f ST < 3 > 

Here AH is a specific type of energy variance [16], 

1/2 



AH = h 



-Oid(O) 



(1) 



2C dd (0) 
and the delta-delta correlation function Cdd(t) is defined [12,16] as 

Cdd (t) = tr ( C -f 3ii ' 2 A a e-^' 2 e ikt / h A b e- iflt ' h ) . (5) 



The generalized delta operator A will be defined below in Eq. (22). 

In applying the QI approximation to the KIEs, it is useful to consider a con- 
tinuous change of the isotope mass. If the two isotopomers are A and B, a real 
parameter A <G [0, 1] can be defined such that 

toj(A) = m Ati (l - A) + m B s\ , 
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where m,A,i and ms,i are the masses of the ith atom in the isotopomcrs A and B, 
respectively. Within the QI approximation (3), the KIE can be expressed as 

KTF - fcQl(0) - Qr(1) x Ag(1) x Cdd(0) x ^W/^dW , fi . 

QI Ml) Qr(0) AiJ(O) C dd (l) C ff (l)/C7 dd (l) ' {> 

where the argument denotes the value of A and for simplicity the time argument 
of the correlation functions has been omitted since it is always t = in the QI 
approximation. Also, for convenience, both numerator and denominator have been 
divided by C dd (A). 

Four types of quantities must be evaluated in order to compute the KIEs from 
Eq. (6): the ratio of the partition functions Q r (l)/Q r (0), the ratio of the delta 
correlation functions C dd (l)/C dd (0), and the energy variance AH (A) and the "ve- 
locity" factor Cff(A)/C dd (A) for A = and 1. The last two quantities are in the 
form of thermodynamic averages (for a given A) and therefore can be directly com- 
puted by Metropolis Monte-Carlo techniques; the relevant estimators have been 
derived in Rcfs. [16, 17]. The most general forms are listed in Rcf. [12]. The first 
two quantities cannot be evaluated directly since they are ratios of quantities for 
two different values of A. 

An elegant solution exists, however. Here is where considering a continuous iso- 
tope change (using a parameter A) becomes useful: instead of computing the ratios 
directly, we use the thermodynamic integration idea [18], applied to the parameter 
A (i.e., to the masses of the isotopes instead to the usual inverse temperature /3). 
We can express the two ratios as an exponential of the integrals of logarithmic 
derivatives, 



Qr(l) _ 

QAO) cxp 

Cdd(l) 

— exp 
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1 dlogQ^A) 

dA dA 
1 dA dlogC dd (A) 



o dA 



(7) 
(8) 



C dd (0) 

Since the logarithmic derivatives can be expanded as 

dlogp(A) _ dp(A)/dA 
dA p(A) ' 

they are normalized quantities (thermodynamic averages) which can be directly 
computed by the Metropolis algorithm. We can compute ratios of both the reactant 
partition functions and the delta-delta correlation functions at A = and 1, by 
computing the values of the corresponding logarithmic derivatives for enough values 
A between and 1, and then by integrating over A and exponentiating, according 
to Eqs. (7) and (8). 

In fact, in a cruder version of the QI method, called the simplest quantum 
instanton (SQI) approximation [11], the ratios of the partition and delta-delta 
correlation functions are all we need, since within that approximation, the kinetic 
isotope effect is just 

KTF - Qril) x Cdd(0) W 

KIESQI QrW) C^T) • (9) 
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The relevant estimators for the logarithmic derivatives have been derived in 
Refs. [12, 13]. In Ref. [12], thermodynamic estimators have been derived that dif- 
ferentiate the kinetic part of the action; in Ref. [13], virial estimators have been 
derived that differentiate the potential part of the action. In both cases, the deriva- 
tions have been done for general systems with TV atoms in d dimensions, even for 
cases with unbound degrees of freedom (such as the center-of-mass motion). In 
the next section, we present a simplified derivation of these estimators for a single 
particle in a one-dimensional external potential. This choice significantly simplifies 
notation, but preserves the main ingredients of the many-dimensional derivation. 

3 Estimators for the logarithmic derivatives of Q r and Cdd 

Below we derive three types of estimators for the logarithmic derivatives of both 
Q r and Cdd- We refer to the three types of estimators as thermodynamic, virial, and 
generalized virial because of their resemblance to corresponding thermodynamic 
[19], virial [20], and generalized virial [21] estimators for the kinetic energy. 

3.1 Partition function 

Let us consider a single particle of mass to in a one-dimensional potential V(r). 
Since we have only one mass to, we do not need to define an additional parameter 
A: we can just take to itself to be the parameter for the thermodynamic integration. 
The PI representation of the partition function is 

M^) P 7 dr<l, '-7 dr( M< rW >). < io > 

p, ({!■<«>}) = exp [-«{rl')))] , (11) 

S— 1 8=1 

Here s = 1, . . . , P, denotes the beads of the discretized paths (s = is identical 
to s = P). In general, we will obtain the estimators for the logarithmic derivative 

directly, by computing the logarithmic derivative — = — — ■ — - of the par- 

dm Q r dm) 

ticular form of the discretized PI. Applying this approach to the PI (10), we obtain 
the thermodynamic estimator 

***-£- -el") , (12) 

dm zto \ am / 



dm 2 / nh 2 [3 2 
4 



E 



r W _ r (s-l) 
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Above, (A({r( s )}) denotes the average over paths weighted with the weight p, 
_JdrW---JdA p )A({r^})p({r^}) 



A 



{-*'})>. 



/ dr(i) • • • / dr( p ) p ({r«}) 
Alternatively, we can define new, mass-scaled coordinates as 

x = m 1/2 r. (13) 

In these new coordinates, the partition function becomes 

M^r/^"-/^"**' (i4 > 



r s=l s=l 



m- 1 / 2 ^ 



Simplest virial estimator for the logarithmic derivative can again be derived by 
direct differentiation of PI (14), 

dlogQ r _ /? /^ dV [m-VV^U = 




dm P \ ^— f dm 

p 



P /* 8V [{m + Am)- 1 / 2 ™ 1 / 2 ^ 8 )] 



Am=0/ 



2P\^ dr^) / 

\s=l / Pr 

Above are three estimators for the logarithmic derivative: the first one suitable if 
the MC simulation is done in mass-scaled coordinates x, the other two for original 
Cartesian coordinates r. The first two suggest evaluation of the derivative numeri- 
cally, by finite differences, which will be in fact, more efficient in many dimensional 
systems than the analytical third expression that requires the knowledge of the 
gradient of the potential. Only for systems with few degrees of freedom and avail- 
able gradient of the potential, the third expression may be preferable. The trick of 
using numerical derivatives with respect to a single parameter was originally used 
by Predescu for computing heat capacities [22] and higher temporal derivatives 
of the flux-flux correlation function [23] where the parameters were the inverse 
temperature and the imaginary time, respectively. 

The simplest virial estimators (15) have one shortcoming compared to the ther- 
modynamic estimators, namely, they only work in bound systems. This can be 
immediately seen by considering a free particle with V(r) = 0. This shortcoming 
can be remedied if the rescaling is done only after subtracting an arbitrarily chosen 
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(but fixed) slice from the remaining P — 1 slices. To be more explicit, let us define 
relative coordinates as 

y (s) = r (s) _ r (p) for s = i ) . . . ) p _ i . 
Since the Jacobian of the transformation is unity, we have 

P/2 



(ssk) J«»-J*r-»J*>".~ 



$ = 



Pm 



" ■' r -U Wy« - y(-l) N 2 



^ +E 



s=2 



2ttS 2 /3 2 

+ 1 ^y(r( p )+y^) + T/(r( p )) 



,(p-i; 



(16) 



Now we define mass-scaled coordinates as 

cc (s) = m 1/2 y (s) - m 1/2 (r W - r (P) ) for s = 1, . . . , P- 1 . 
In these coordinates, the partition function becomes 

P/2 



Vr 



r = (— P — \ m V2 f dx (D . . . f dx (P-D f dr W e -/J* , 



$ = 



\2irh 2 (3J 
P 



yOlVvLW.^-iAVf^-i)^ 



2tt?i 2 /3 2 

+ 1 ^T/(r( p )+m- 1 /2 x ( s )) + T/(r( p )) 



(17) 



(18) 



The generalized virial estimator for the logarithmic derivative follows by differen- 
tiating the PI expression (18), 



dlogQ r _ 1 /3 /A dV [r( p ) + m-^W] 



dm 2m P \-^^ 9m 

x s— 1 



1 £ / p ^ [r( p ) + (m + Am)- 1 / 2 ™ 1 / 2 (r(*) - r< p ))] 
2^ ~ P\2^ 



8=1 
P 



<9Am 



Am=0/ p 



1 g /f^ OV(rW) (a) r(P) . 

9-m OP\Z-^ Q-f»1 v 



2m ' 2P\^ 9r( s ) 

X S — 1 



(19) 

Since we have chosen the slice s = P arbitrarily, we can do the same for any slice 
s, derive a corresponding estimator, and then take an average of these estimators. 



The result is 
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dlogQ 
dm 



= - + -(E 



f3 /Pdv(A s i) 



dr (s 



(r 



W 



(20) 



s=l 



.(«) 



and in general, we can replace r^ ' in all three forms (19) of the estimator by r c . 

Since the number of slices P appears explicitly only in the denominator of the 
generalized virial estimator (19) or (20), the statistical error should be independent 
of P for a fixed number of Monte-Carlo samples. On the other hand, P appears 
explicitly in the numerator of the thermodynamic estimator (12), so the error is 
expected to grow with P. This will be confirmed in the numerical example in 
Section 4. 

3.2 Delta— delta correlation function 

The derivation for Cdd is similar. However, due to the constraints to the two 
dividing surfaces, a new term appears in the estimator. The PI representation of 
Cdd is 



C, 



.id 



P % ({r (s) } 



(S) /^•••/d^y({^}), (2D 

& (r<°>)] A [$, (A p ^) 



The generalized delta function A is defined [16] as 



A[£(r)] 



drt 



5[£(r)] 



(22) 



For numerical purposes, it is convenient to replace the strict delta function by a 
Gaussian approximation [16], 



A 



e r 
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(23) 



r («) + r (»+i) 



A[£(r)] = 



2P 

tt/i 2 /? 



1/2 



cxp 



2Pto 

' h 2 (3 



[frac£(r)d r €(r)]' 



We can define an effective action $ e fj = $ + K nstr , which includes the constraint 
potential 



V a 



2Pm 

h 2 fi 2 

7 



Mir). 



(24) 
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The logarithmic derivative will have one extra term due to this constraint potential, 



d log C d d d log Q r 



dm 



dm 



dV c , 



dm 



(25) 



where — : (p r — > p*) denotes that the corresponding estimator for Q r given in 

dm 
Eq. (12), (15), or (19) should be used except that the sampling is done according to 

weight p* instead of p r . Using the PI representation (21) of Cdd m Cartesian coor- 
dinates, we find the additional term from Eq. (25) to the thermodynamic estimator 
(12) to be 



dV cc 



'IP 



£(r) 



H 2 f3 2 



cUM 



dm h 2 (3 2 

Rescaling coordinates according to Eq. (13) gives a constraint potential 

2 
7Pm. flmr'-^x) 

V„ 



(26) 



£ (m 1 / 2 x) 
d r £ t (m~ l l 2 x) 



The additional term from Eq. (25) to the simple virial estimator (15) becomes 

2 



dV cx 



2P 



dm 



H 2 f3 2 dAm 



(rn + Am) 



1/2 



S V rcscj 

Cr r e S cS V fCSC ) 



(27) 



Ar 



Finally, if we rescalc coordinates according to Eq. (17), or better, as 

x (s) =m 1/2 (r (s) -r c ) for s = 1, . . . , P , 
we obtain the same estimator as (27), only the rescaled coordinate is defined as 

1/2 



-(«) 



Ar 



( r W_ r c). 



(28) 



The generalization to more-dimensional systems is fairly straightforward. Only 
in the case of the simple virial estimator (15) or (27), care must be taken to ac- 
count for the unbound (free) degrees of freedom by appropriately rescaling the 
corresponding volume. For instance, for bimolccular reactions, the potential in the 
reactant region is independent of the center-of-mass-coordinate and the relative 
coordinate of the two molecules. For details, see Rcf. [13]. 



4 Numerical results 

In Ref. [12], the QI procedure for evaluating KIEs was successfully tested on sev- 
eral problems of increasing complexity: the one-dimensional Eckart barrier and the 
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isotopic variants of both the collinear and the three-dimensional hydrogen-exchange 
reaction H+H 2 -► H 2 +H. The results for the KIE = fc(H+H 2 )/fc(D+D 2 ) as a func- 
tion of the inverse temperature \/T for both collinear and three-dimensional ver- 
sions of the reaction are also shown here in Fig. 1 . The figure compares the exact 
quantum-mechanical result [12, 24] with the results of the QI, SQI, and TST ap- 
proximations. The three-dimensional version also shows the result of the canonical 
variational TST with semiclassical tunneling correction (CVT) [25] . In general, the 
results of the QI approximation are very good: the error is smaller than 10% for 
temperatures 250 to 600 K. For lower temperatures, the larger error is due to us- 
ing a single dividing surface: results can be improved by considering two separate 
dividing surfaces. At high temperatures, the error is due to classical recrossing 
and cannot be corrected within the QI approximation. For further details of the 
calculation see Rcf. [12]. 




2 3 

1000 / T(K) 



2 3 

1000 / T(K) 



Fig. 1. Kinetic isotope effect A;(H+H2)/fc(D+D2) for the hydrogen exchange reaction: 
(a) the collinear version, (b) the three-dimensional version. 



Three types of estimators for the logarithmic derivative of Q r , given in Eqs. (12) , 
(15), and (20), are compared in Fig. 2. This calculation is for the collinear version of 
the KIE = £;(H+H 2 )/fc(D+D 2 ) at 300 K. The calculation was done with 10 walkers 
and a fixed number 10 5 Monte-Carlo moves for all P. The left part of the figure 
shows the convergence of the partition- function ratio Q r (l)/Qr(0) as a function of 
the number of slices P. The ratio was computed via the thermodynamic integration 

(7) in which the three different estimators (12), (15), and (20) for 



dA 



were 



used. The right panel shows the P-dcpcndcnce of the relative error of Q r (l)/Q r (0). 
As expected, for large P, the error of the generalized virial estimator (20) is almost 
independent of P. On the other hand the error of the thermodynamic estimator (12) 
grows with P. Even for small P, the generalized virial estimator is superior. Finally, 
we can see the importance of subtracting the centroid motion before rescaling in Eq. 
(17) by comparing errors of the simple (15) and generalized (20) virial estimators. 
The difference is due to the fact that we have two free degrees of freedom in the 
reactant region of the collinear bimolecular reaction. Similar conclusions (not shown 
here) can be obtained for the ratio Cdd(l)/Cdd(0), except that the error of the 
generalized virial estimator has a weak dependence on P arising from the additional 
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Fig. 2. Comparison of the thermodynamic, virial, and generalized virial estimators for 

the logarithmic derivative of Q r . Left: ratio of the partition functions as a function of P, 

right: its relative error as a function of P. 



term due to the constraint to the dividing surfaces. 

5 Conclusion 

Judging from the numerical results in the previous section, the QI approach 
for computing kinetic isotope effects is very promising. The procedure is quite 
general: multi-dimensional estimators for all relevant quantities are presented in 
Refs. [12,13]. Although the path-integral approach has been chosen because of 
its favorable scaling with the number of degrees of freedom, the computations for 
many-dimensional systems are still difficult. 

One obstacle is the difficulty of efficient sampling of a many-dimensional config- 
uration space. With proper estimators, we may decrease statistical and systematic 
discretization errors, but it is difficult to avoid systematic errors due to long corre- 
lations. For this reason, it may be efficient to use a different number of imaginary 
time slices [26] for different degrees of freedom, which is a generalization of more 
crude mixed quantum-classical methods. 

Another obstacle to obtaining a good match between theory and experiment 
is the potential energy surface for the reaction. While more accurate ab initio 
potentials are computationally very expensive, the much faster molecular-mechanics 
force fields are often too crude. In Ref. [13], in which the QI method is used to 
compute the KIE for the isomcrization of cis-pentadiene, two approaches are taken: 
in one, an empirical valence bond (EVB) potential is formed from the equilibrium 
potentials for reactants and products; in the other, a more accurate but also a more- 
expensive semi-empirical potential is used. Because of the computational expense 
already for this system with 39 degrees of freedom, it appear that due to their 
better accuracy such semi-empirical potentials will be the potentials of choice for 
intermediate-size systems, and the EVB potentials based on molecular-mechanical 
force fields the potentials of choice for truly many-dimensional systems. 
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